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Abstract 

The density matrix renormalization group is a variational approxi- 
mation method that maximizes the partition function — or minimize 
the ground state energy — of quantum lattice systems. The variational 
relation is expressed as Z = Tr p > TV (lp), where p is the density sub- 
matrix of the system, and 1 is a projection operator. In this report we 
apply the variational relation to two-dimensional (2D) classical lattice 
models, where the density submatrix p is obtained as a product of the 
corner transfer matrices. The obtained renormalization group method 
for 2D classical lattice model, the corner transfer matrix renormaliza- 
tion group method, is applied to the q = 2 ~ 5 Potts models. With 
the help of the finite size scaling, critical exponents (q = 2, 3) and the 
latent heat (q = 5) are precisely obtained. 
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1 Introduction 



The basic procedure in the renormalization group (RG) is to keep relevant 
information of a physical system, and neglect (or integrate out) irrelevant 
one. jy, || |3| The density matrix renormalization group (DMRG) introduced 
by White [[§] greatly enhances the applicability of the numerical RG, because 
the method automatically keeps a fixed numbers (= m) of the relevant basis; 
DMRG present the best approximation within the limited numerical resource 
that we can use. The DMRG has been applied to a number of one-dimensional 
(ID) quantum lattice systems, such as the spin chain, || |6| ladder, |7], |j Bethe 



lattice system, || strongly correlated electron systems, [jTT], [12|, [13], mod- 
els in momentum space. Not only the numerical superiority, but also the 
formulation of DMRG have attracted theoretical interests. Ostlund and Rom- 
mer jnj have shown that DMRG is a variational method, where the ground 
state is expressed as a product of 3-index tensors. [|T^, [T7| Martm-Delgado 
and Sierra have investigated the analytic formulation of DMRG, and have for- 
mulated the correlated block RG. [|18| |19| Recently White have refined the 
finite system algorithm of DMRG, and extend the applicability of DMRG to 
2D quantum systems. p0[ Quite recently Xiang have reported DMRG study 



of ID quantum system at finite temperature, 21] using the quantum transfer 



matrix formulation |22] and DMRG applied to the transfer matrix. [23 



Another RG approach has been development for 2D classical lattice models: 
Baxter's method of the corner transfer matrix. (CTM) |24| The method is 
a generalization of Kramers- Wannier approximation, |25], and therefore 
Baxter's method is based on a variational principle for the partition function. 
It should be noted that Baxter's variational relation is in principle the same as 



the variational relation in DMRG. |27| The purpose of this report is to explain 
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how the concept of DMRG is applied to 2D classical lattice models. We start 
from a short review of the variational relation in DMRG in the next section. 
It is worth looking at a practical use of the RG method as the 2D photo 



image compression. [28, B^J A photo image in our computer is normally com- 
pressed before it is stored, in order to decrease the file size. The compression 
algorithm is related to the block RG method, fl], H since a small region — a 
pixel — in a 2D picture has strong correlation with its environment. (Com- 
puter scientists may insist that their findings about the photo image compres- 
sion are efficient for RG formulation in physics.) At present, compression of 



movies — TV pictures — are in progress in the world of computation; [3C , 31 
the algorithm may be a good reference for the RG study of 3D classical systems 
and 2D quantum systems. 

The development of RG, DMRG, Baxter's CTM method, and the photo 
image compression is summarized in Fig.l. Originally these methods are pro- 
posed independently, however, now it is apparent that their background is 
in common. It is, to approximate a system (or the objects) within a limited 
number of freedom. 

2 Variational Principle in DMRG 

We start from a short review of the variational principle in DMRG. We consider 
the antiferromagnetic S — 1/2 Heisenberg spin chain as an example of ID 
quantum systems. The spin Hamiltonian is 

F = j£>-S m , (1) 

i 

where S; represents the quantum spin at i-site, and the parameter J is positive. 
The Hamiltonian H is real-symmetric, and so is the density matrix p = e~^ H . 
(In the following discussion, p does not always have to be real-symmetric, but 



should be positive definite.) 

We consider an open spin chain, (Fig. 2) which consists of the left half [L] 
(= the local system) and the right half [R] (= the reserver). The terms 'local 
system' and 'reserver' are rather formal, since [R] is not always longer than 
[L\. (In Fig. 2 both [L] and [R] has the same size.) The Hilbert space of the 
whole system is spanned by the real-space basis 1 1 ) | r) , where 1 1 ) and | r) 
corresponds to the spin configuration for [L] and [R], respectively. The matrix 
element of p is given by 

Pir,Vr'^(l\(r\e^ H \r')\l'). (2) 

In the context of DMRG, what is called the density matrix is actually the 
density submatrix (DSM) 

/=ziop«<<n = £io(£/^v)<n (3) 

IV IV V r ' 

that contains the information only about the local system [L\. The trace of p L 
is equal to the partition function. 

The relevant state selection — renormalization — in DMRG is performed 
through the diagonalization of the DSM 

OVQ = diag{A 1 ,A 2 ,...}, (4) 

where Ai > A 2 > • • • > are eigenvalues in decreasing order, Q = (qi, q 2 , . . .) 
are the set of the corresponding right eigenvectors, O = (o!,o 2 , . . .) are that 
of the left eigenvectors. The matrices O and Q satisfies the dual orthogonal 
relation QO T = 1. When p L is real- symmetric, Q is equal to O, and both of 
them are orthogonal matrices. The first m column vectors q a (1 < a < m) 
in Q represent the RG transformation from the original basis {| l) \ r)} to the 
renormalized basis {| a)}- The irrelevant states are thrown away. 0] 

3 



We can check the validity of the above basis state selection by observing 
the inequality 

m 

z = Ea,>5>. (5) 
i 1=1 

The quantity Z = YaLi \ is the approximate partition function, which is 
smaller than Z by J2i> m M- It has been known that if the system has a fi- 
nite excitation gap, the eigenvalue Aj decays exponentially with respect to i, 
and therefore Z is a good approximation for Z when m is sufficiently large. 
The approximate partition function is equal to the trace of the m-dimensional 
diagonal matrix 

p L = diag{Ai, A 2 , . . . , A m } = 6 T p L Q, (6) 

where Q is the rectangular matrix (q x , q 2 , . . . , q m ), and O is (o 1? o 2 , . . . , o m ). 
We can regard the matrix operation of T and Q to p L as the RG transfor- 
mation (or the block spin transformation), and p L as the renormalized DSM. 

Substituting Eq.(6) into Eq.(5), we obtain a variational relation in the 
matrix form 

Z = Tr p L > Trp L = Tr (QO T p L ), (7) 

where the matrix product 1 = QO T has the property of the projection operator 
(QO T ) 2 = QO T . The projection operator QO T is optimal in the sense that 
it gives maximum of Z = Tr(lp) under the constraint l 2 = 1 and Trl = m. 
In other word, the DMRG minimize the free energy of the system within the 
restricted degree of freedom. At the zero temperature ((3 — > oo), DMRG 
minimize the total energy. 

The optimal projection operator 1 for the local system [L] is dependent 
on the size of the reserver [R]. However, the dependence is not conspicuous 
when [R] is sufficiently large. The infinite system DMRG algorithm uses the 
insensitivity of 1 against the reserver size. The finite system DMRG algorithm 
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is more accurate than the infinite algorithm, because the former correctly takes 
into account of the reserver-size dependence. 

3 From Quantum system to Classical system 

In order to apply the variational relation in DMRG to 2D classical lattice sys- 
tem, we define the DSM for 2D systems. We use the fact that the density 
matrix e~^ H of a d- dimensional quantum system can be expressed as a par- 
tition function of a <i+l-dimensional classical system with special boundary 
conditions; the relation is known as the Trotter-Suzuki formula. |32|, |33| 
The density matrix of the Heisenberg chain in Fig. 2 is approximated as 

e -P« = (e-^Y ~ (e-$ H *e-$ a °y (8) 

where Ha = J2i &2i • S 2 j+i and H B = J2i S 2 ,,+i • S 2i+2 are the partition of H. 
The matrix element pi r ^ T > in Eq.(2) is approximated by the Boltzmann weight 
of the chessboard model, whose boundary spin configurations are fixed to I, r, I' 
and r' . Figure 3 shows an example when M — 4. When we consider the 
partition function Z, the boundaries Z, r, V and r' plays the role of the tabs for 
sticking; Z = Tr e~^ H is approximately equal to the partition function of the 
cylindrical system shown in Fig.4(a), which is constructed by attaching / and 
r in Fig. 3 to I' and r', respectively. In the same way, the approximate DSM 
p\ v = J2 r Pir,i'r is obtained by attaching r to r' as shown in Fig. 4(b). 

The expression of p L in Fig. 4(b) is a typical example of the DSM for 2D 
classical lattice models. The p L corresponds to the cylindrical system with a 
cut L, where the spin configurations around the cut are fixed to I and I'. We 
generalize this example. Suppose we have a finite size lattice system [A] shown 
in Fig. 5. We then consider an arbitrary line or curve L on [A], and cut [A] 
along L to derive a new system [A']. (The curve L is a kind of string on the 
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1+1 space time.) The derived system [A'] has new boundaries around the cut 
L, where the boundary spin configurations are represented by the labels I and 
V. The Boltzmann weight of [A'], which we write pf v , is the DSM of [A]. 

We have defined p L for 2D classical lattice models. Once we obtain p L for 
a given system [A], we can perform DMRG along the variational treatment 
discussed in the previous section, where the group of spins I and V on [A] are 
transformed into m-state effective spins. In the next section we present an 
appropriate choice of [A] and L for a typical 2D classical model. 

4 Construction of the Density Matrix via CTM 

The Trotter-Suzuki decomposed Heisenberg spin chain in Eq. (8) is an anisotropic 
2D lattice model, and the application of DMRG on this system is rather 
complicated. [^TJ For the tutorial purpose we consider a simpler 2D system, the 
symmetric 16-vertex model. p4| , |33| The model includes the Ising model |35| 



as its special case. For simplicity, we assume that the model is ferromagnetic. 
The 16-vertex model is defined by the Boltzmann weight 

W abcd = W bcda = W cdab = W dabc (9) 

on each vertex (= lattice point) of the simple square lattice, where the spin 
variables a,b,c, and d take either + (up) or — (down). In the following we 
consider a square system [A], whose linear dimension is 2N or 2N + 1. We 
impose fixed boundary condition on [A]; we introduce boundary weights 

Pabc = W abc+ (10) 

and 

C ab = W ab++ (11) 
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in order to fix the boundary spins to +. (Fig. 6) The partition function of [A] is 
expressed by these weights. For example, the partition function of the system 
with 2N + 1 = 3 is expressed as 

Z = W kheb P abc C cd P de fC '^ g P ghi C^Pj kl C lal (12) 

ab...l 

where the position of spin indices a-d is shown in Fig. 7. 

In order to generalize Eq.(12) to arbitrary system size, we introduce a half- 
row transfer matrix, (HRTM) which is a generalization of the boundary weight 
P in Eq.(lO). The HRTM of length N is defined by the recursion relation 

P&bc = ^ <u d cjy b P&'d c' i (13) 

d 

where the label N is the number of vertices in HRTM, P^ bc is equal to P abc in 
Eq.(lO), and a represents a group of spins on a row 

a = (ai, a 2 , • • • , On_i, a N ), (14) 

which is related to a' as a = (a', ajv); the same for c = (c', cjy). Figure 8(a) 
shows an example when N = 3. We occasionally drop the vector indices of P^ bc 
and write it simply as P b N ; in that case we think of P b as a 2 Ar -dimensional 
matrix (P b N ) ac - The HRTM is also conventionally called vertex operator. [^] 
The density submatrix p L of the square system [A] is expressed as a product 
of corner transfer matrices. (CTMs) 0] The CTM of size N is defined as 

C^b — ^2 [ ^2^efb N a N P&'e c' ^b' / d' ) ^c'd^ ( ^) 

c'd' V ef 7 

where we have used the index rule a = (a', ajv), b = (b',6jv), etc. Figure 
8(b) shows the example when N = 2. The smallest CTM C* b is equal to the 
boundary weight C ab in Eq.(ll). The square system [A] is then constructed 
by attaching four CTMs. Figure 9 shows the system [A'] of the size 2N = 6, 
whose Boltzmann weight 
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corresponds to the DSM of the square system [A] of the size 2N = 6. Such a 
construction of the DSM was first introduced by Baxter more than 30 years 
ago in his variational method. 



Now we can apply the RG procedure to 2D classical lattice models using p L 
defined in Eq.(16); following the RG procedure in Sec. 2, we combine Baxter's 
method and DMRG. For the brevity, we call our new RG method corner 



transfer matrix renormalization group (CTMRG) in the following. [27 

5 CTM Renormalization Group 

As was discussed in Sec. 2, the heart of DMRG is the diagonalization of p L . 
For the symmetric 16-vertex model, where p L is expressed as (C N ) 4 , both p L 
and C N have the common eigenvectors. We therefore diagonalize CTM 

T C N Q = diag{^,^ 2 ,...} (17) 

instead of p L , where we assume the decreasing order \ui\ > \u 2 \, ■ ■ ■ > 0. 
The block-spin transformation is performed by the rectangular matrices O = 
(oi, o 2 , . . . , o m ) and Q = (qi, q 2 , . . . , q m ), where o a and q Q are the left and 
the right eigenvectors of C , respectively. Here after we use greek letters 
for indices that runs from 1 to m. (Actually, Q is equal to O since C N is 
symmetric.) The renormalized CTM 

C N = T C N Q = diag{^, u 2 , . . . , u m } (18) 

is related to p L = T p L Q via p L = (C N ) 4 . The renormalized HRTM is 
obtained in the same way 

P a N = T P a N Q, (19) 
where the tensor elements of are P£L- 
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At this point we remember that P N and C N are defined through the recur- 
sion relations Eq.(13) and Eq.(15), respectively. The relations are also valid for 
P N and C N . For C N , its area is extended by attaching two HRTMs (Fig.10) 

C(a } a)l/3,b) = ^2 ^efba ^aeS Ppf5 u S, (20) 
efS 

where C,£yip h s is a partially renormalized CTM of linear size N + 1, and 
the pairs of indices (a, a) and (/3, b) represent the row and the column matrix 
indices for C N+l . The length of HRTM is simultaneously increased by putting 
a new vertex at the end point 

-^(a,a)\( 7 ,c) = W adcb Pad~{ ■ (21) 
d 

The extended CTM in Eq.(20) is not diagonal. As was done in Eqs. (17- 
18), we diagonalize C + to obtain the new renormalized CTM C N+1 with 
the matrix dimension m. It is now obvious that we can indefinitely repeat 
the renormalization process in Eqs. (18-19) and the system size extension in 
Eqs. (20-21). In this way, we obtain P N and C N for arbitrary large N starting 
from P 1 and C 1 . [|7| 

Approximate thermodynamic functions of [A] can be obtained using C . 
The free energy —ksT In Z is estimated from the approximate partition func- 
tion Z = Tr (C N ) 4 = J^i^t- K is also possible to obtain local quantities, such 
as spin polarization and multi-spin correlation functions, using a combination 
of P N and C . For example, the local energy is estimated as 



EY TW piV s-iN tdN r<N tdN r<N v»N 'r*N\ 
v N abed ^ abed 11 {^a U r b U r c U r d U ) ^2) 

Tabe d w abcd i:T{p^c^p b N c^p^c^p^c^y 

where the position of the indices a-d are shown in Fig. 11, and X abc d is the 
local energy operator hx(W a b c d)- Since the CTMRG extends the area of the 
system [A] from the center, local quantities at the center can be calculated 
most precisely. 



The largest matrix element of C rapidly grows with respect to N. We 
should therefore normalize CTM 

d " diag{l,^...,^U(7" (23) 



L0\ V UJ\ U>i 

in realistic numerical calculations. We should also normalize P^ in the same 
manner. Apart from the critical point, the normalized C converges to its 
thermodynamic limit C°° exponentially with respect to N. At criticality the 
convergence is relatively slow; it is observed that the decay rate at criticality 



is controlled by the critical exponents. 37 



We have imposed fixed boundary conditions on [A]. Since the boundary 
condition is totally determined by the boundary weight P\ in Eq.(lO) and C 1 
in Eq.(ll), we can choose other boundary conditions by modifying P^ and C . 
For example, the free boundary condition is imposed by the boundary weights 

Pabc 

W abc+ + W abc - (24) 

and 

C ab = W ab++ + W ab+ . + W ab „ + + W ab — (25) 

The CTMRG method presented above can be applied to a wide class of 2D 
classical lattice models, such as the g-state Potts model, IRF model, etc. We 
have to be careful to anisotropic lattice models whose p L consist of four differ- 
ent CTM's. (In Baxter's textbook the DSM are expressed as p = ABCD. 0) 
In such a case, we should not diagonalize each CTM independently, but we 
should diagonalize p L to perform RG transformation. We should also be care- 
ful to antiferromagnetic models, because we have to prepare several sets of P N 
and C N according to the alternating spin order. 
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6 Numerical Result 



The symmetric 16-vertex model includes the square lattice Ising model as its 
parameter limit. We first look at the critical phenomena of the Ising model 
using CTMRG. Figure 12 shows the calculated local energy 

E(2N) = Tr {aa'p L )/Z (26) 

where aa' is a pair of neighboring spin at the center of the square cluster [A] 
of size 2N. The data shown by the black dots are E(oo), that are obtained 
when m = 98. The data deviate from the exact solution at most 10 -7 . At 
the critical temperature T c , we estimate E c (oo) by observing its convergence 
with respect to N. The inset of Fig. 11 shows the 1/N dependence of E C (N). 
A simple 1/N fitting gives E c (oo) = 0.707148, which is close to the exact one 
l/y/2 = 0.707107.... 

The 1/N dependence of E C (N) is actually related to the scale invariance 
at criticality. The finite size scaling (FSS) theory |38|, ^ predicts that E C (N) 
obeys the scaling form 

E C (N) - E c (oo) ~ N^- d , (27) 

where v is the correlation-length exponent and d — 2 is the spatial dimension- 
ality; in our case E C (N) is proportional to iV -1 because v of the Ising model 
is equal to unity. Similarly, the local order parameter 

M(N) = Tr (ap L )/Z (28) 

at the center of [A] obeys 

M C (N) ~ AH d - 2 +")/ 2 (29) 

with the anomalous dimension of the spin r\. Figure 13 shows the iV dependence 
of calculated magnetization M C (N) at T c when m = 148. As it is expected 
from r) = 1/4, the calculated M C (N) is proportional to N^ 1 ^ 8 . 
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Table 1: Critical exponents v and r/ of the q=2, 3 Potts model estimated from 
the numerical data when m = 200. Theoretical values are shown inside the 
parenthesis. 



q 


v (Exact) 


f] (Exact) 


2 
3 


1.0006 (1.0000) 
0.8321 (0.8333) 


0.2501 (0.2500) 
0.2654 (0.2667) 



The estimation of the critical exponents via FSS is valid for a wide class 
of 2D classical model. As examples, we apply CTMRG and FSS to the q=2, 3 
Potts models, [50, 51] which can be treated as a symmetric g 4 -vertex model. [52, 



4341 Table I shows the estimated exponents v and r\ from the calculated data 



when m = 200 and 100 < N < 1000. For comparison, the theoretically 



determined exponents [4T| are shown in the parenthesis. The error in the 
calculated exponents are less than 0.2%. 

Numerical analysis of the case q = 4 is in progress. In this case, one has to 



take care of the logarithmic corrections in the FSS analysis. [[|3] 



We finally show the N dependence of the local energy E(N) of the q=5 
Potts model at the transition temperature; we measure E(N) from the average 
E* = (E + + E-)/2, where E + and E_ are the energy of the disordered and 
the ordered phase, respectively, at the transition temperature. (Fig. 14) Since 
the transition is first order, E(N) — E* does not obey the scaling form in 
Eq.(27). In this case, it is possible to estimate the latent heat through the 
scaling analyses with respect to both N and m. The calculated latent heat is 



0.0256, where the exact value is 0.0265. |41] The CTMRG is expected to be 



efficient for the analysis of the weak first order transition. J44 

7 Conclusion and discussion 



We have reviewed the variational principle in DMRG, and explained its ap- 
plication to 2D classical lattice systems. We create the DSM according to 
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Baxter's construction, and perform the RG transformation using DMRG al- 
gorithm. As trial calculations, we perform FSS analyses on the g-state Potts 
models at the transition temperature. It is concluded that the CTMRG is 
efficient for the determination of the critical exponents or the latent heat. 

It is possible to feed-back the CTMRG algorithm to the zero-temperature 
ID quantum system, and to obtain a rapid infinite system DMRG algorithm; 



the product wave function (PWF) RG. [§6| The PWFRG is closely related to 
the improved finite system DMRG algorithm (DMRG++ in Fig.l) recently 
proposed by White. f20fl 

Finally, we discuss the applicability of DMRG or CTMRG to 3D classi- 
cal lattice models. As we have employed a square system for 2D models, we 
should consider a cubic system for 3D models. We divide it into 8 subcubics, 
say corner transfer tensor, (CTT) and construct the DSM as their tensor prod- 
ucts. We then diagonalize the DSM, and renormalize CTT. Apart from the 
renormalized CTM in two dimension, renormalized CTT is not diagonal. [33] 



There is, however, a computational problem that prevents us from practical 
use of the RG method for 3D system; at present the numerical calculation is 
too heavy even for a small m. 

There are plenty of subject that we can clarify in the field of the numerical 
RG. 
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Figure Captions 

Fig. 1. Historical overview of DMRG and related fields. 

Fig. 2. The S — 1/2 Heisenberg spin chain with open boundary conditions. We 
divide it into the local system [L] and the reserver [R]. 

Fig. 3. Trotter-Suzuki decomposition of the density matrix in Eq.(8). The label 
/, r, /' and r' denote the boundary spin configurations. 

Fig. 4. The density matrix in the Trotter- Suzuki formula Eq.(8): (a) partition 
function Trp. (b) density sub matrix p\ v . 

Fig. 5. We derive the DSM of the finite size system [A] by cutting it along the 
curve L. Compared to [A], the system [A'] has additional boundaries I 
and 

Fig. 6. Boltzmann weights of the symmetric 16-vertex model; P and C are the 
boundary weights in Eq.(lO) and Eq.(ll), respectively. 

Fig. 7. A square system of the linear dimension 2N + 1 = 3, whose partition 
function is given by Eq.(12). 

fig. 8. Recursive definition of (a) P N in Eq.(13) and (b) C N in Eq.(15). The 
shown examples are P 3 and C 2 . 

fig. 9. The density submatrix p L of the square system [A] is expressed as a 
product offourCTMs. (Eq.(16).) 

Fig. 10 Area extension of the renormalized CTM in Eq.(20). 

Fig. 11. A square system of size 2N + 1 is expressed as a product among C N , 
Pa% cd) , and W. 
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Fig. 12. Local Energy E(oo) of the Ising Model at the center of the square system. 

The inset shows the size dependence of E(N) at the critical temperature 
T 

Fig. 13 Local Magnetization of the Ising Model at T c . 
Fig. 14 Local Energy of the q=5 Potts model. 
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